﻿#pragma execution_character_set("utf-8")
#include <stdexcept>
#include "TheTask.h"
#include "gdal_priv.h"
#include <vector>
using std::vector;
using std::runtime_error;

#define THETASK_NODATA -9999

TheTask::TheTask(string infileB1, int band1Based1,
	string infileB3, int band3Based1,
	string outfile, double coefA, double coefB, double scale, double coefC)
{
	_infileB1= infileB1;
	_b1Based1= band1Based1;
	_infileB3= infileB3;
	_b3Based1= band3Based1;
	_outfile= outfile;
	_coefA= coefA;
	_coefB= coefB;
	_coefC = coefC;
	_scale = scale;
}
void TheTask::run()
{

	string error;
	GDALDataset* indsB1 = 0;
	GDALDataset* indsB3 = 0;
	GDALDataset* outds = 0;
	try {
		indsB1 = (GDALDataset*)GDALOpen(_infileB1.c_str(), GA_ReadOnly);
		if (indsB1 == 0) {
			throw runtime_error("打开输入B1文件失败。");
		}
		indsB3 = (GDALDataset*)GDALOpen(_infileB3.c_str(), GA_ReadOnly);
		if (indsB3 == 0) {
			throw runtime_error("打开输入B3文件失败。");
		}
		int xsize = indsB1->GetRasterXSize();
		int ysize = indsB1->GetRasterYSize();
		int xsize1 = indsB3->GetRasterXSize();
		int ysize1 = indsB3->GetRasterYSize();


		if (xsize == 0 || ysize == 0 || xsize!=xsize1 || ysize!=ysize1 ) {
			throw runtime_error("xsize, ysize 其中一个为0或者B1波段与B3波段图像尺寸不一致。");
		}
		GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
		outds = driver->Create(_outfile.c_str(),xsize,ysize,1, GDT_Float32  ,0);

		double trans[6];
		indsB1->GetGeoTransform(trans);
		outds->SetGeoTransform(trans);
		outds->SetProjection(indsB1->GetProjectionRef());

		doOne(indsB1, indsB3, outds, _b1Based1, _b3Based1, xsize, ysize, _coefA, _coefB,_scale,_coefC);
		outds->GetRasterBand(1)->SetNoDataValue(THETASK_NODATA);
		if (indsB1) GDALClose(indsB1); indsB1 = 0;
		if (indsB3) GDALClose(indsB3); indsB3 = 0;
		if (outds) GDALClose(outds); outds = 0;
		this->notifyDone();
		emit onSaveTaskDone();
	}
	catch (runtime_error& ex) {
		if (indsB1) GDALClose(indsB1); indsB1 = 0;
		if (indsB3) GDALClose(indsB3); indsB3 = 0;
		if (outds) GDALClose(outds); outds = 0;
		error = ex.what() ;
		this->notifyError(error);
	}
}

void TheTask::doOne(GDALDataset* indsB1, GDALDataset* indsB3, GDALDataset* outds, int b1B1, int b3B1,
	const int xsize, const int ysize,
	double cA, double cB, double scale, double cC)
{
	//processing
	int per0 = -1;
	double perScale = 100.0 / ysize;
	vector<float> rowdata1(xsize);
	vector<float> rowdata2(xsize);
	vector<float> rowdataOut(xsize);
	for (int iy = 0; iy < ysize; ++iy) {
		indsB1->GetRasterBand(b1B1)->RasterIO(GF_Read, 0, iy, xsize, 1, rowdata1.data(), xsize, 1, GDT_Float32, 0, 0, 0);
		indsB3->GetRasterBand(b3B1)->RasterIO(GF_Read, 0, iy, xsize, 1, rowdata2.data(), xsize, 1, GDT_Float32, 0, 0, 0);
		for (int ix = 0; ix < xsize; ++ix) {
			double r1 = rowdata1[ix]*scale;
			double r2 = rowdata2[ix]*scale;
			if ( 0.0<r1 && r1<1.0 && 0.0<r2 && r2 <1.0 ) {
				float c = cA * (r1 - r2) + cB;
				if (c > 1.0) {
					c = cC / sqrtf(c);
				}
				else {
					c = cC;
				}
				rowdataOut[ix] = c;
			}
			else {
				rowdataOut[ix] = THETASK_NODATA;
			}
		}
		outds->GetRasterBand(1)->RasterIO(GF_Write, 0, iy, xsize, 1, rowdataOut.data(), xsize, 1, GDT_Float32, 0, 0, 0);
		int per1 = (iy+1) * perScale;
		if (per1 != per0) {
			per0 = per1;
			this->notifyProgress(per0, 100);
		}
	}
}